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Abstract. The identification and characterization of massive galaxies over a wide redshift range allow us to place stringent 
constraints on the cosmic history of galaxy mass assembly and on current models of galaxy formation and evolution. This 
work explores the existence of high redshift massive galaxies unveiled with Spitzer+IRAC, but missed by conventional selec- 
tion techniques based on optical and near-infrared observations. To this end, we use multi-wavelength imaging data available 
for the GOODS-South field (130 arcmin^), and select a flux-limited sample from the IRAC 3.6 /jm image to 53 6 ^ 1.8/jJy 
(m(AB)<23.26). In order to identify the most extreme objects and to complement previously published selections in this field, 
we confine our study to the galaxies undetected by the optical HST-l-ACS imaging and close to the detection limit of the A'-band 
image (K > 23.5 AB). Our selection unveiled 20 galaxies on which we performed a detailed analysis. For each galaxy, we built 
a Spectral Energy Distribution (SED) based on optical-to-8/jm photometry. The SEDs were then used to estimate the photo- 
metric redshifts and to derive the main galaxies physical properties. Further constraints were also obtained from the available 
X-ray and 24yum data. The majority of the sample (14 out of 20) sources show degenerate/bimodal solutions for the photometric 
redshifts. These can either be heavily dust-enshrouded (Ay ~ 2 - 4) starbursts at 2 < z < 3 with bolometric luminosities 
LiR > 10'"Lo, or massive post-starburst galaxies in the redshift interval 4 < z < 9 with stellar masses of ~ 10" Mq. The 
remaining six galaxies present a less ambiguous photometric redshift: with the exception of one low-z dusty source, these latter 
objects favour a low-extinction solution, with four of them showing best-fit photo-z solutions at z ~ 4. One galaxy, ID-6, the 
only source in our sample with both an X-ray and a 24 /jm detection, might be an extremely massive object at z ~ 8 detected 
during a post-starburst phase with concomitant QSO activity responsible for the 24 yum and X-ray emissions (although a lower-z 
solution is not excluded). Our investigation of Spitzer-selected galaxies with very red SEDs and completely undetected in the 
optical reveals a potential population of massive galaxies at z > 4 which appear to include significant AGN emissions. These 
sources may be the oldest stellar systems at z ~ 4, given that the estimated ages are close to the age of the Universe at that red- 
shift. We found that these, previously unrecognized, optically obscured objects might provide an important contribution to the 
massive-end (M > 10" Mq) of the high-z stellar mass function and they would almost double it. Our suggested evidence in these 
mature high-z galaxies of the widespread presence of hidden AGNs may have important implications for galaxy formation, due 
to their feedback effects on the surrounding ISM. 



1. INTRODUCTION 

Early attempts to estimate the luminosity and mass functions 
of galaxies have revealed a surprisingly low rate of evolution 
of the stellar mass between the present epoch and redshifts ~1, 
and particularly so for the most massive systems (Fontana et 
al. 2004; Bundy et al. 2005; Franceschini et al. 2006). Deep 
extensive near-IR surveys have recently confirmed the presence 
of a numerous population of already massive galaxies at z ~ 
2-4 (e.g., Franx et al. 2003; van Dokkum et al. 2003; Cimatti 
et al. 2004; Daddi et al. 2004; Yan et al. 2004; Le Fevre et al. 
2005; Papovich et al. 2006). 
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However, very little is currently known about the exis- 
tence, number density and properties of massive galaxies be- 
yond z ~ 4. Mobasher et al. (2005) have recently reported 
the discovery of a very massive (M > 10" Mq) galaxy in the 
GOODS-South field at z ^ 6.5, although a lower-z solution has 
been suggested by Dunlop et al. (2006). Using the Early Data 
Release (EDR) by the Ultra Deep Survey (UDS) component of 
the UKIRT Infrared Deep Sky Survey (UKIDSS), McLure et 
al. (2006) have found nine Lyman-break galaxy candidates at 
z ~ 5 -6 with stellar masses larger than 5 x lO'^M©. The forma- 
tion epochs of such massive galaxies, and of their stellar pop- 
ulation contents, cover a wide range of redshifts, whose upper 
boundary extends well into the re-ionization epoch (Panagia et 
al. 2005). 
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Fig. 1 Plots of the K - [3.6] versus z - [3.6] colours for our reference GOODS/CDFS IRAC-selected sample with 5 3.6 > 1.8 yuJy 
(small black dots in the top panel). The objects from our high-redshift sub-sample correspond to the big red circles with z-band 
limits. Data are compared with the colours of five spectral templates: a 10 Gyr old and passive elliptical (solid line), a star-forming 
(dotted line) and a post-starburst galaxy (dashed), a single Simple Stellar Population (SSP) with an age of 1 Gyr (unextinguished, 
three-dots-dashed line) and a younger extinguished SSP of 100 Myr (Ay - 2.4, dot-dashed line). The predicted colours are shown 
as a function of redshift starting from z - with increasing steps of 0.5 (the z = values fall typically in the bottom-left part of 
the panel). The bottom panel zooms onto the colours of our candidate very-high-redshift objects undetected in the z-band down 
to z = 28.15 (Icr limit), and an additional sample of objects fainter than z - 26.95 (3cr limit). 



A complementary information is offered by very high- 
redshift quasars (Fan et al. 2003; Fan et al. 2004) located in 
highly metal-enriched interstellar environments (e.g. Dietrich 
et al. 2003; Freudling et al. 2003). These cases are a clear in- 
dication of already advanced evolutionary stages of the host 
galaxies at such high redshifts (e.g. Maiolino et al. 2006 and 
references therein), particularly considering that the identified 
atomic species are subject to delayed supernova enrichment, 
requiring activity at very high redshifts. Further information 
comes also from the detection of thermal dust continuum and 
molecular emissions in these quasars including the farthest 
known QSO at z = 6.42, SDSS Jl 148-1-5152, indicating the 
presence of large amounts of dust (M > 10** Mq) and molec- 
ular gas (M > lO''* Mq) in their circumnuclear environments 
(Omont et al. 1996; Omont et al. 2003; Bertoldi et al. 2003; 
Robson et al. 2004). The formation of dust grains requires con- 



densation processes in addition to the stellar activity required 
to produce the basic elements (C, Si), and in normal condi- 
tions occur in the expanding envelopes of evolved stars (AGB, 
late giants, with evolutionary timescales of > 10^ yrs). If the 
dust heating came from starburst activity, as it has often been 
suggested, that would imply an enormous star formation rate 
(SFR) of several hundred to several thousand solar masses per 
year in the objects. 

2. OBSERVATIONS AND SAMPLE SELECTION 

Finally, the nuclear black-hole mass estimates of the highest-z 
SDSS quasars range from several times 10^ to several times 
lO'' Mq (e.g. Fan 2006), such that even with continuous 
Eddington-limited accretion they should have started to form 
at z > 10. Although direct dynamical estimates may suggest a 
kind of break of the relation between the black hole and the host 
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Fig. 2 Multiwavelength identification of the twenty 3.6 /im-selected objects. For each source we report a postage of 5"x5" in the 
ACS z-band, of 8"x8" in the ISAAC and IRAC bands and of 17"xl7" in the 24 fim MIPS band. North is up, East is at left. 



galaxy mass at the highest redshifts (Walter et al. 2004; Peng 
et al. 2006), it is hard to imagine that such supermassive black 
holes do not reside in some kind of forming massive galactic 
bulge. 

In conclusion, both directly detected high-z galaxies and 
studies of the metal enriched circumnuclear media in quasars 
imply enhanced star formation activity at very high redshifts 
(z > 6) in some specific cosmic sites. Establishing how fre- 
quent the latter are in representative cosmic volumes would set 
important constraints on galaxy formation models, particularly 
considering the fast decline of the dark matter halo mass func- 
tion with increasing redshifts. Indeed, the number densities and 
mass functions of host galaxies at the highest redshifts are far 
from settled at the moment, while lunited information is avail- 
able only on the most luminous quasars (e.g. from SDSS & 
2DF). 

High redshift galaxies are currently selected using a vari- 
ety of techniques. A very successful approach is based on the 
detection of the spectral break in the UV continuum blueward 
of the Lya due to the intervening Lya-forest. By construction, 
this technique is biased towards star-forming galaxies with rest- 
frame UV fluxes not strongly reddened by dust extinction (e.g. 
Steidel et al., 2003; Bouwens et al. 2006). This selection tech- 
nique is biased against selecting high-redshift galaxies with red 
spectra, either due to dust extinction or to the presence of ma- 



ture stellar populations. Thus, other approaches have been re- 
cently used to select high-z galaxies in a way less affected by 
biases. A successful example is represented by near-IR surveys 
with or without additional color or photometric redshift selec- 
tions (e.g. Cimatti et al. 2002; Franx et al. 2003; Abraham et al. 
2004; Daddi et al. 2004), or pure flux-limited optically-selected 
samples with no color cuts (e.g. Le Fevre et al. 2005), or sub- 
millimeter/millimeter selection of dusty galaxies (Small et al. 
2002, Dannerbauer et al. 2004). 

The main question is then: are there high redshift galax- 
ies that are missed by the current optical and near-IR selection 
techniques? 

The advent of Spitzer Space Telescope opened a new pos- 
sibility in this respect as it allows to select samples in a spec- 
tral region (3-8jum) that was not accessible from the ground. 
The selection of samples at these wavelengths becomes par- 
ticularly important for the specific case of high redshift mas- 
sive galaxies because the 3-8//m selection allows to sample 
the rest-frame near-IR of the high-z galaxy Spectral Energy 
Distributions (SEDs) and is therefore sensitive to their stellar 
mass rather than to their star formation activity, and also much 
less affected by dust extinction effects. 

In this work, we explored if there are high redshift mas- 
sive galaxies missed by the conventional selection techniques, 
but that can be unveiled with Spitzer. For this purpose, we 
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Fig. 3 Postage stamps showing the result of the stacking procedure of all ACS bands (B, V, /, z) at the positions of the twenty 
near-IR dark sources. No optical signal is detected. The postage has a size of 5" x 5". 
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Fig. 4 K - [3.6] colour versus K magnitude plot for our candidate high-redshift galaxies (filled circles) and for galactic stars of 
different spectral types, as indicated. K magnitudes are calculated by setting stars at a 30 kpc distance. 



used the multi-wavelength deep imaging data available for the 
GOODS-South field (130 arcmin^, GiavaHsco et al. 2004) and 
we searched for extreme galaxies in a way complementary to 
the other selections applied in previous surveys in the same 
field. 

The paper is organized as follows. In Section 2 we present 
the multi-wavelength photometric data-set and our sample se- 
lection criteria. In Section 3 we discuss the model fitting to 
the observed SEDs, while in Section 4 we analyze the infrared 
colours of our sample. Section 5 illustrates the X-ray properties 
of the sample sources. In Section 6 we report the detection of a 
candidate massive galaxy at z ~ 8. Section 7 presents a general 
discussion of the statistical properties of our sample. Finally, 
in Section 8 we summarize the main results of this paper We 
adopt for the cosmological parameters Q„,= 0.3 and Qa = 0.7, 
and H() = 70 Km/s/Mpc. All magnitudes are given in the AB 
system. 

With the aim of studying the multiwavelength photometry 
of extragalactic sources measured by different instruments, we 
need to measure the bulk of the emission from each object in 



each photometric band. Only with such kind of approach an 
SED can be considered reliable. In our work, we performed 
aperture photometry in each band. 

2.1. Deep S pitzer Near- and Mid-IR Photometric 
Imaging 

As part of the GOODS project, the Spitzer Space Telescope has 
recently surveyed the CDFS field in the IR between 3.6 and 8.0 
jum using IRAC and in the range 24-160 //m using MIPS. The 
fully reduced data were publicly released by the GOODS team 
and are available through the World Wide WetQ. 

In this paper we exploit a galaxy catalogue that we have de- 
rived from the 3.6 //m IRAC public raw data and recently used 
to derive the luminosity and mass function by morphological 
type (Franceschini et al. 2006). In order to obtain the most ac- 
curate SEDs, we have recomputed the fluxes of each source 
by performing aperture photometry in the four IRAC bands at 

' http://data.spitzer.caltech.edu/popular/goods 
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Fig. 5 J-K colour versus /T- [4.5] colour for our candidate high-redshift galaxies (filled circles). The 7-band magnitudes adopted 
in this plot correspond to the 2cr value (J - 26.29). The dashed region indicates the expected colours for Galactic brown dwarfs, 
as recently suggested by Mannucci et al. (2006). 



the positions originally detected in the IRAC 3.6 /zm channel. 
Assuming that essentially all the sample sources are seen as 
point-like by the IRAC ~ 2 arcsec FWHM PSF imager, we 
computed with SExtractor (Bertin & Arnouts, 1996) the fluxes 
within a 3.8 arcsec diameter aperture. This choice is supported 
by an accurate analysis performed by the SWIRE tearrB They 
constructed color-magnitude diagrams for various types of ob- 
jects, in particular main-sequence stars. It was then found that 
the scatter in these diagrams is minimized through the use of 
a 3.8 diameter aperture, and corresponds to roughly twice the 
beamwidth. To obtain total fluxes, we then applied the correc- 
tion factors indicated by the SWIRE teanJl. We have indepen- 
dently verified that the IRAC/SWIRE aperture corrections are 
consistent with those derived by fitting the radial brightness 
profiles of stars in the GOODS fields. In the case of extended 
sources, we used Kron-like magnitudes (AUTO_MAG output 
parameter in SExtractor). Our sample turned out to be ~ 60% 
complete above 1 /iJy (m3.6=23.9), ~ 75% complete above 2 
//Jy (m3.6=23.15), ~ 90% at 5 /iJy (m3.6=22.15), and more than 
~ 95% above 10 juJy (m3.6=21.4). 

The MIPS public dataset includes calibrated maps and a 
catalogue of 24 jiva sources with flux densities ^24 > 80yuJy. 
The photometry is based on a PSF fitting algorithm, where the 
SExtractor positions of the IRAC sources are used as input to 
the MIPS source extraction process. The MIPS 24 //m PSF was 
generated from isolated sources in the image, and renormalized 
based on the aperture corrections published in the MIPS Data 
Handbook (v2.1, section 3.7.5, table 3.12). 



^ http://data.spitzer.caltech.edu/popular/swire/20050603_enhanced_v 
^ See Note 2 



To extend the 24 yum sample to fainter fluxes, we have run 
an independent PSF fitting algorithm that we akeady success- 
fully applied in the GOODS -xFLS/ENl Science Verification 
field (Rodighiero et al. 2006). By these means we have ex- 
tended the 24 fim. sample down to 5 24 > 20juJy. 

2.2. Near-IR Ground-based Imaging 

As part of GOODS, near-infrared imaging observations of the 
CDFS have been carried out in the J, H and bands, using 
the ISAAC instrument mounted on the ESO VLT. We made 
use of the publicly available J, H and Ks imaging (version 1 .0, 
releasee!] by the ESO/GOODS team in April 2004). This data 
release includes 21 fully reduced VLT/ISAAC fields in 7, H 
and K,, covering 130 arcmin^ of the GOODS/CDFS region. It 
also includes mosaics of the coadded tiles as single FITS files in 
J and Ks bands, as well as the corresponding weight maps. To 
provide a homogeneous photometric calibration across the en- 
tire field, all images are rescaled to the same zero point (26.0). 
The final mosaics have a pixel scale of 0.15". 

We measured the J, H and K band magnitudes with 
SExtractor at the positions of the 3.6 IRAC selected 
sources through circular apertures with diameters of 2 arcsec. 
Mobasher et al. (2005) found that the photometric curve of 
growth converges at this aperture, which represents the best 
compromise between convergence of the total flux and the ef- 
fects of systematic uncertainty in the background subtraction. 

For undetected ISAAC sources, we have initially computed 
the upper limits to the flux by measuring the signal in 400 ran- 
dom sky positions. We used 2-arcsec aperture diameters and 

1/ 

* http://www.eso.org/ science/goods/releases/20040430/ 
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Fig. 6 The IRAC colour [8.0] - [4.5] is shown against the colour [5.8] - [3.6]. The dashed line marks the region dominated by 
AGNs in the redshift range 1 < z < 3 according to Lacy et al. (2005). The near-IR dark objects of our sample are plotted as 
red filled circles. For comparison, we report a sample of SCUBA radio-selected sources detected at Spitzer wavelengths (open 
diamonds, Frayer et al. 2004) and a sample of Chandra sources detected in the SWIRE survey by Franceschini et al. (2005). We 
plot the evolutionary tracks for a reddened starburst (Arp220 - green solid lines), a Seyfert 1 galaxy (dashed cyan lines) and a 
passive elliptical (dot-dashed red lines). The predicted colours are shown as a function of redshift with increasing step size of 
0.5, starting from z - 0.5 (as marked in the plot for each template). For comparison, we also report the colour distribution (dotted 
points) of the GOODS/MUSIC sample (Grazian et al. 2006) with 5(8jum) > ljuJy. 



calculate the value of the standard deviation from the distri- 
bution of the measured aperture fluxes. We obtained Icr val- 
ues of 26.12, 25.82 and 25.12 magnitudes for the /, H and 
K band, respectively. However, given that the depth of the 
ISAAC imaging varies significantly across the CDFS field, we 
preferred to compute independent fluxes for each undetected 
source. Following the approach of Dunlop et al.(2006), we per- 
fomed our own manual photometry through a 2-arcsec diame- 
ter aperture at the IRAC-selected positions. 

2.3. ACS/HST Optical Imaging 

The core of the GOODS project was the acquisition and data 
reduction of high-resolution HST/ACS imaging obtained as an 
HST Treasury Program (Giavalisco et al. 2004). The GOODS 
ACS/HST Treasury Program has surveyed two separate fields 
(the CDFS and the Hubble Deep Field North) with four broad- 
band filters: F435W(B), F606W(V), F775W(i) and F850LP(z). 
In August 2003 the GOODS team released version 1 .0 of the 



reduced, stacked and mosaiced images for all the data ac- 
quired over the five epochs of observation. To improve the point 
spread function (PSF) sampling, the original images, which had 
a scale of 0.05 arcsec / pixel, have been drizzled on to images 
with a scale of 0.03 arcsec / pixel. 



The dataset is complemented by the ACS/HST catalogues 
released by the HST/GOODS team in October 2004, contain- 
ing the photometry in B, V, i and z bands. The source extraction 
and the photometric measurements have been performed by the 
GOODS team running a modified version of SExtractor on the 
version 1.0 images. We have considered aperture magnitudes 
using a 1 arcsec diameter For undetected IRAC sources, pho- 
tometry in the ACS B, V, i and z bands was measured as upper 
limits with the same procedure adopted for the ISAAC imag- 
ing. We used 1 arcsec diameter apertures and adopted the 3cr 
value as an upper limit to the flux of undetected sources. 
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Fig. 7 The colour [24]-[8] is shown as a function of the colour [8]-[2.2]. Our sources are plotted as big filled circles. We also 
report evolutionary tracks for starbursts galaxies (Arp220 - green solid lines), AGNs (Seyfert 1 - dashed cyan lines, Seyfert 2 
- dashed yellow line), combined quasar/ultraluminous infrared galaxy (ULIRG) sources like the type-1 QSO Mrk 231 (orange 
solid Unes) and the type-2 quasar/ULIRG Superantennae South (IRAS 19254s; Berta et al. 2003, magenta lines), and a passive 
elliptical (dot-dashed red lines). The predicted colours are shown as a function of redshift with increasing step size of 0.5, starting 
from z - 0.5 (as marked in the plot for each template). 



Table 1 Photometric properties of the sample objects 



ID 


RA 


DEC 


b 


V 


1 




z 


J 


H 


K 


[3.6] 


[4.5] 


[5.8] 


[8.0] 


[24] 




J2000 


J2000 


mag 


mag 


mag 


mag 


mag 


mag 


mag 


mag 


mag 


mag 


mag 


mag 


1 


53.12440 


-27.88268 


-27.61 


-27.72 


-27 


15 


-26.94 


-25.41 


-25.65 


24.49 


22.72 


22.35 


22.13 


21.62 


18.69 


2 


53.22620 


-27.85910 


-27.61 


-27.72 


-27 


15 


-26.94 


-26.19 


-24.97 


23.95 


22.58 


22.32 


21.92 


22.20 


19.03 


3 


53.06607 


-27.83178 


-27.61 


-27.72 


-27 


15 


-26.94 


-25.78 


-25.17 


23.67 


22.22 


22.06 


21.94 


22.08 


19.27 


4 


53.19752 


-27.81387 


-27.61 


-27.72 


-27 


15 


-26.94 


-26.86 


-25.06 


-24.86 


23.17 


22.85 


22.83 


21.82 


18.78 


5 


53.16146 


-27.81118 


-27.61 


-27.72 


-27 


15 


-26.94 


-25.45 


24.94 


24.15 


22.11 


21.81 


21.60 


21.37 


19.63 


6 


53.15831 


-27.73355 


-27.61 


-27.72 


-27 


15 


-26.94 


-25.65 


24.68 


24.42 


22.22 


21.65 


21.32 


20.85 


19.27 


7 


53.15552 


-27.71871 


-27.61 


-27.72 


-27 


15 


-26.94 


-25.37 


-25.06 


-24.30 


23.26 


23.02 


22.71 


22.93 


20.70 


8 


53.19850 


-27.92742 


-27.61 


-27.72 


-27 


15 


-26.94 


-26.27 




23.92 


23.23 


22.86 


22.64 


21.87 


00.00 


9 


53.13935 


-27.89070 


-27.61 


-27.72 


-27 


15 


-26.94 


-26.65 


-25.92 


25.11 


22.73 


22.13 


21.52 


21.46 


19.69 


10 


53.14963 


-27.87682 


-27.61 


-27.72 


-27 


15 


-26.94 


-25.37 


-25.30 


-25.12 


22.97 


22.64 


22.42 


22.33 


19.72 


11 


53.14664 


-27.87103 


-27.61 


-27.72 


-27 


15 


-26.94 


-25.37 


-26.30 


24.09 


23.03 


22.47 


21.99 


21.44 


00.00 


12 


53.12938 


-27.87172 


-27.61 


-27.72 


-27 


15 


-26.94 


-25.39 


-25.20 


-24.95 


23.00 


22.57 


22.26 


21.85 


00.00 


13 


53.13809 


-27.86780 


-27.61 


-27.72 


-27 


15 


-26.94 


-26.79 


-26.02 


24.15 


23.14 


22.82 


22.68 


22.64 


00.00 


14 


53.20188 


-27.84421 


-27.61 


-27.72 


-27 


15 


-26.94 


-25.55 


-24.88 


24.37 


23.10 


22.76 


22.51 


22.35 


19.40 


15 


53.03842 


-27.82529 


-27.61 


-27.72 


-27 


15 


-26.94 


-26.05 


-25.08 


24.26 


22.89 


22.56 


22.35 


22.39 


00.00 


16 


53.18194 


-27.81414 


-27.61 


-27.72 


-27 


15 


-26.94 


-26.01 


-24.82 


-24.08 


23.02 


22.72 


22.10 


20.84 


17.78 


17 


53.12701 


-27.80455 


-27.61 


-27.72 


-27 


15 


-26.94 


-25.54 


-25.38 


24.04 


22.74 


22.55 


22.39 


22.46 


00.00 


18 


53.18047 


-27.77972 


-27.61 


-27.72 


-27 


15 


-26.94 


-25.54 


-24.91 


23.70 


22.25 


22.10 


21.80 


21.86 


00.00 


19 


53.11983 


-27.74309 


-27.61 


-27.72 


-27 


15 


-26.94 


-26.25 


-26.02 


24.67 


22.28 


21.80 


21.32 


20.50 


00.00 


20 


53.12762 


-27.70684 


-27.61 


-27.72 


-27 


15 


-26.94 


-25.37 


-24.82 


23.82 


22.62 


22.41 


22.33 


21.92 


00.00 



* Source #8 is out of the H band image. 
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Fig. 8 Left panel: confidence levels for the photometric redshift z and extinction Ay derived from ?lx^ analysis are shown models 
for each source, for the BC03 model. Dotted, dashed, dot-dashed and three-dot-dashed curves respectively mark the 68%, 90%, 
95% and 99.99% confidence levels of the;^'^ statistics. For each source, the red shaded area shows regions of the parameter space 
which appear to be disfavoured by our analysis in Sect. 14. II taking into account the 24 yum constraint. Right panel: the value of 
the is plotted as a function of redshift. Different curves show the result of using different extinction ranges in the SED fitting 
procedure with Hyperz (solid curve: Ay < 6, red dashed curve: Ay < 10). 
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2.4. Sample Selection 

Our main aim is to search for galaxies possibly missed by pre- 
vious surveys and to find new candidates of high redshift mas- 
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sive systems. In a recent paper, Dunlop et al. (2006) made use to investigate the high-z galaxy population in the GOODS-S 
of a /T-band selected sample brighter than K - 23.5 AB mag 
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field. They conclude that there is no convincing evidence for In order to search for extreme galaxies "dark" at optical 
any galaxies with M > 3 10" at z > 4. wavelengths and ultrafaint in the near-lR, and to be comple- 
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mentary with the selection applied by Dunlop et al. (2006) in 
the same field, we selected our sample starting from the Spitzer 
IRAC 3.6 jjm image, and applied the following selection crite- 
ria: 

(1) flux-limited selection from the IRAC 3.6 /im image 
down to ^3.6 > 1.8/zJy (m(AB)<23.26). 

(2) objects undetected in any of the optical HST-nACS im- 
ages at the lo" level (i.e. z > 28.14) 

(3) objects close to the detection Umit of the VLTh-ISAAC 
/T-band deep image (K > 23.5 AB, a complementary criterion 
to that of Dunlop et al. 2006, including objects with K < 23.5). 

In order to avoid potential photometric uncertainties from 
the blending due to the large IRAC PSF, we excluded from 
our sample a few blended sources. We found twenty sources 
satisfying the above conditions. This final sample was further 
investigated on ISAAC images: the twenty sources turned out 
to be all undetected in the /-band (at least at the ISAAC sen- 
sitivity), while only two candidates showed a faint counterpart 
in the //-band image. Three sources are undetected even in the 
K-hand ISAAC image. 

To illustrate the eff'ects of our selection criteria, we report in 
Figure[T]the K - [3.6] versus z - [3.6] colours for our complete 
GOODS/CDFS sample with 5 3.6 > 1.8 yuJy (top panel). For 
comparison, we also show here the colour-colour plots of five 
spectral template SEDs: a 10-Gyr old passive elliptical, a star- 
forming and a post-starburst galaxjH, a solar-metallicity SSP 
with an age of 1 Gyr (unextinguished) and a younger extin- 
guished SSP of 100 Myr (Ay - 2.4). The predicted colours are 
shown as a function of redshift starting from z = with in- 
creasing steps of 0.5 (the z - value falls in the bottom left 
part of the panel). The bottom panel zooms onto the colours 
of our candidate very-high-redshift objects undetected in the 
z-band down toz = 28.15 (Icr limit), and, for comparison, an 
additiona sample of objects fainter than z - 26.95 (3cr limit). 
The meaning of the symbols is detailed in the caption of Figure 

m 

We note that ~ 40 sources, falling in the colour-colour re- 
gion populated by our final sample, escaped our selection crite- 
ria being all brighter in the K- and z-bands with respect to our 
imposed constraints (K > 23.5 and z > 28.14). 

A summary of the multiwavelength identifications of the 
3.6 fim objects is presented in Figure |2] for each sources we 
report here postage stamps of 5" x 5" size in the ACS z-band, 
8" X 8" in the ISAAC and IRAC bands and of 17" x 17" in the 
24 ijm MIPS band. 

In order to look for any possible faint optical detection, we 
stacked together the images of all sample sources in the four 
ACS bands. The result of this stacking procedure is presented 
in Figure [3] showing the complete absence of optical signal at 
the position of the selected sources. 



' The synthetic spectra are taken from the set of template used in 
Fritz et al. (2007, submitted). The first is taken as representative of a 
galaxy during a post-starburst phase with a second main episode of 
star formation -forming ~ 10% of the total stellar mass- at ~ 10** 
years while the second is build with a main-impulsive burst at ~ 13 
Gyr and a continuous star formation which is truncated at 2 • 10' years. 



The photometric data used to construct the spectral energy 
distributions discussed in the next Sects, are presented in Table 
1 . Negative values correspond to upper limits. We adopt a com- 
mon value of 10% (15%) of the measured fluxes as photometric 
eiTors for the IRAC (MIPS) bands, in order to reflect the sys- 
tematic uncertainties of the instruments The main contributions 
to these uncertainties are due to the colour-dependence in the 
flat field and to the absolute calibration (see for example Lacy 
et al., 2005, and the IRAC and MIPS Data Handbook). For the 
ISAAC K- and //-band fluxes, we adopt as photometric error a 
standard value of 0.3 magnitudes (as derived from the standard 
deviation statistics from our set of Monte Carlo simulations dis- 
cussed in Franceschini et al., 2006). 

The observational SEDs for all our sample sources are plot- 
ted in Figs. |9] and below. 

2.5. Contamination by Galactic Stars 

We have first considered the possibility that some sources 
within our sample could be misidentified Galactic cool stars. 

We verified that ultra-cool galactic stars (like M and L 
dwarfs) show colours (e.g. K - [3.6], see Figure |4| that are 
in general much bluer than those of our sources. Moreover, 
the peak of the stellar emissions of M dwarfs falls shortward 
(A ~ I - 2yum) than observed in the SEDs of our sampl^ 

Another possible contaminant is represented by evolved 
dusty stars, like AGBs (in particular carbon stars). The spec- 
tra of these objects can in principle reproduce the SEDs of our 
sample in the 0.4-8.0 yum interval range (see again Figure |4]l. 
However, the detection of several AGB stars in an area of only 
130 armin^ at the magnitude limit of K < 23.5 is much in ex- 
cess of what is e.g. observed in the LMC (Cioni et al., 2006). In 
addition, such objects are so bright that they should be moved 
far away from our Galaxy (at least -15 Mpc) in order to match 
the observed fluxes of our sample objects. 

Young stellar objects (YSOs) may also be considered. 
Although their spectra resemble those of our sources, they are 
usually associated with extended molecolar clouds in Galactic 
star-forming regions that are obviously absent in the GOODS- 
South field. 

An interesting possibility is the potential contamination by 
Galactic brown dwarfs. Mannucci et al. (2006) have recently 
suggested that a pair of faint z-dropout sources (z - J > 0.9) in 
the CDFS are compatible with the expected colours of brown 
dwarfs. We checked this hypothesis in Figure |5] where we re- 
port the J - K colour versus the K - [4.5] colour for our candi- 
date high-redshift galaxies. The 7-band magnitudes appearing 
in this plot correspond to our 2cr value (J - 26.29). The dashed 
region indicates the expected colours for galactic brown dwarfs 
(adapted from Mannucci et al., 2006, that used the stellar mod- 
els of Allard et al., 2001). Clearly, all sources in our sample 
present J - K colours much redder than expected for brown 
dwarfs. 



^ We made use of the colour predictions from the stel- 
lar model of Jarrett et al. (1994). We also used the 
spectral templates of the ISO standard stars available at 
ihttp://www.iso.vilspa.esa.es/users/expl_lib/ISO/wwwcal/, 
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Fig. 9 Best-fit models for the sample sources with a bimodal solution. In the left panel we report the low-redshift solution. The 
observed SED of each source (open diamond) is shown together with the corresponding low-redshift best-fit solution (solid red 
line) indicated by Figure [8] and Table 3 for the stellar component (up to 8 jum). The IR starburst template that better reproduces 
the S {24 jjin) / S i^i^m) flux ratio is reported as a dashed cyan line. The IR spectra have been normalized to match the 24 /im 
measurements. In the right panel we shows the observed SEDs together with the corresponding best-fit high-z solutions (dashed 
red lines, sol. II in Table 3) In this case, the IR part of the spectra have been reproduced with the spectral template of a dusty torus 
representing the emission of a type-2 AGN (green dotted lines, Fritz et al. 2006 model). The IR spectra have been normalized to 
match the 24 /im measurements. The solid red lines correspond to the sum of the different galaxy components. 
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Figure |9] (continued) . 



3. INFRARED COLOURS 

3.1. Near-IR to IRAC Colours 

To start constraining the nature of our selected galaxies, we 
have first compared their mid-lR properties with recently sug- 
gested colour-colour diagnostic plots. While the composite 
spectra of the stellar populations in normal galaxies produce 
SEDs peaking at approximately 1 .6 yum, the U V to mid-IR con- 



tinua of AGNs are dominated by power-law emission. Based 
on this. Lacy et al. (2004) used IRAC colours from the Spitzer 
First Look Survey to identify AGNs. A region in the [4.5]-[8.0] 
versus [3.6] - [5.8] colour plot where AGN are expected to lie 
is shown as a dashed-line contour in Figure |6l 

The twenty near-IR dark objects of our sample are plot- 
ted as red filled circles. For comparison, we also report the 
colour distribution (dotted points) of the GOODS/MUSIC sam- 
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Fig. 10 Best-fit models for the sample source with a single pho- 
tometric solution. The observed SED of each source (open di- 
amond) is shown together with the corresponding best-fit so- 
lution (dashed red line) indicated by Figure |8] for the stellar 
component (up to 8 fim). The IR part of the spectra have been 
reproduced with the spectral template of a dusty torus repre- 
senting the emission of a type-2 AGN (green dotted lines, Fritz 
et al. 2006 model). The IR spectra have been normalized to 
match the 24 jjm measurements. The solid red lines correspond 
to the sum of the different galaxy components. 



pie (Grazian et al. 2006) with SiS/jm) > lyuJy. The colour 
tracks for a dusty starburst (Arp220 - green solid lines), a 



Seyfert 1 galaxy (dashed cyan lines) and a passive elliptical 
(dot-dashed red lines) are also shown as a function of red- 
shift, with increasing step size of 0.5, starting from z - 0.5 
(as marked in the plot for each templates). All sources fall 
within the AGN area. However, this is not necessarily an in- 
dication of an AGN dominance in our sample. Indeed, Fig. 
|6] also shows the colour distribution of a sample of SCUBA 
radio-selected sources detected at Spitzer wavelengths (open 
diamonds, Frayer et al. 2004), and the SWIRE/Chandra sam- 
ple selected in the ENl field by Franceschini et al. (2005, as- 
terisks). 

Clearly, both starforming galaxies, ellipticals and AGNs at 
z > 2 can reproduce the IRAC colours of our sample, showing 
a limited diagnostic power by the test. 



3.2. Near-IR -to- mid-IR colours 



It is remarkable to note that half of our 3.6 fim selected sources 
in our sample (11 out of 20) reveal a significant 24 //m excess. 
Given the MIPS and IRAC limiting fluxes, this is inconsistent 
with a purely passive galaxy at any redshift, and is a clear indi- 
cation for some activity taking place in the objects. The nature 
of this activity (star formation, AGN or the concomitance of 
both processes) will be subject of investigation in the present 
work. 

In Figure|2]the colour [8] -[24] is shown as a function of the 
colour A'- [8] (see also Webb et al. 2006), where our sources are 
plotted as big filled circles. We also report here colour tracks for 
starbursts galaxies (Arp220 - green solid lines), AGNs (Seyfert 

1 - dashed cyan lines, Seyfert 2 - dashed yellow line), a com- 
bined quasar/ultraluminous infrared galaxy (ULIRG) sources 
like the type-1 QSO Mrk 231 (orange solid lines), the type- 

2 quasar/ULIRG Superantennae (IRAS 19254s; Berta et al. 
2003, magenta lines), and a passive elUptical (dot-dashed red 
lines). These templates are derived from the spectral library re- 
ported in PoUetta et al. 2006 (see also Franceschini et al. 2005). 
The predicted colours are shown as a function of redshift with 
increasing step size of 0.5, starting from z = 0.5 (as marked in 
the plot for each template). 

The extremely red infrared colours of our sample are dif- 
ficult to explain. The K - [8] colour on the X-axis might be 
consistent with the SED of a starburst galaxy, like Arp220, or 
alternatively with a passively evolving galaxy, both at z > 3. 
Lower-redshift solutions require additional extinction (Ay ~ 2) 
to that of the Arp220 spectrum. The very red [8] - [24] ratio 
may be consistent alternatively with a moderate-redshift (z < 3) 
dusty starburst (Arp220), or with a dusty quasar. 

The combined set of colours are reproduced only by high 
redshift (z > 2.5) ULIRGs with concomitant QSO activity (the 
type-2 quasar Superantennae, more marginally by the infrared- 
luminous type-1 QSO Mrk231). This colour-colour plot of Fig. 
|2]will be considered in later applications. 
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4. SPECTRAL ENERGY DISTRIBUTION ANALYSIS 

4.1. Synthetic Spectral Models and Photometric 
Redshift Estimates 

We made use of the Hyperz code (Bolzonella et al. 2000) to 
estimate the photometric redshifts of each source. The whole 
broad-band photometric dataset available to us was exploited, 
with the exclusion of the 24 flux (the stellar population syn- 
thesis model template spectra do not include dust emission). 
The fitting procedure was based on a maximum-likelihood al- 
gorithm and the quality of the fit is investigated by means of a 

statistics. The code computes the for a given number of 
templates, which differ for star formation histories, metallici- 
ties and ages, and finds the best-fitting template among them. 

We used two stellar population synthesis models to fit the 
observed SEDs: Bruzual & Chariot (Bruzual & Chariot 2003, 
hereafter BC03) and the Maraston (2005, hereafter MA05) 
models. In both cases, we assumed exponentially decreasing 
star formation rates (SFR) parameterized by a time-scale t. We 
have considered the following set of values for t: 0. 1, 0.3, 1, 2, 
3, 5, 15, 30 Gyr. Moreover, we included a case with constant 
SFR, and a single burst model corresponding to an individual 
simple stellar population (SSP). For the BC03 model we used 
a Chabrier Initial Mass Function (IMF). The MA05 evolution- 
ary tracks available to us have been generated with a Salpeter 
IMF. To convert mass estimates obtained with Salpeter IMF 
to Chabrier IMF we used a constant value of 0.23 dex (which 
corresponds to correct the masses computed with Salpeter IMF 
downwards by a factor 1.7). This value has been obtained by 
comparing the stellar masses estimated with BC03 models built 
with different IMFs for an observed spectroscopic sample of 
galaxies. In this way we obtained the mass difference expected 
for galaxies with the same colours. 

The extinction parameter is allowed to span the widest 
conceivable range of values < Ay < 10 (with a step of 
d{Av) =0.1). We assumed the extinction law by Calzetti et 
al. (2000), and the metallicity was set to the solar value, in 
order to minimize the number of free parameters at play. We 
also accounted for Lyman-series absorption due to HI clouds in 
the intergalactic medium, following the prescription of Madau 
(1995). Even if the reddening is formally allowed to reach ex- 
treme values, we remind that values of Ay exceeding ~5-6 
magnitudes have been observed only in the inner part (i.e. the 
central 1-2 kpc) of local luminous IR galaxies (e.g. Mayya et 
al. 2004, Poggianti et al. 2001). In comparison, at high red- 
shift (2 < z < 3) the dusty submillimeter galaxies have typ- 
ical average extinction around Ay ~ 2.4 (Small et al. 2004, 
Kjiudsen et al. 2005). Another class potentially including heav- 
ily obscured sources is that of Distant Red Galaxies (DRGs) 
and dusty Extremely Red Objects (EROs). However, even in 
this case the reddening has typical values around Ay <~3.0 at 
2 < z < 3 (Cimatti et al. 2003, Moustakas et al. 2004, Stern 
et al. 2006, Papovich et al. 2006). An Hyper Extremely Red 
Object (HERO) has been proposed to lie at z=2.4 with Av=4.5 
mag (Im et al. 2002). However, no convincing evidence for 
galaxies with a global value of Ay exceeding 5-6 mag has been 
reported until know. 



4.2. Additional constraints from 24^m flux 

With the attempt of reducing the uncertainties in the extinction 
estimate, we have considered the MIPS 24 //m flux in addition 
to the 0.4-8 /im photometric data. The dust emission spectrum 
longwards of few jim (restframe) has been modelled by assum- 
ing the observed IR SED of Arp220, and then compared to the 
flux detected in the MIPS band. The Arp220 IR template cor- 
responds to the spectrum of an highly absorbed ultra-luminous 
IR starburst, a conservative choice for our analysis. 

The predicted IR emission for each one of our objects is 
calculated as the difference between the unextinguished and 
the extinguished optical emission, assuming that all the flux ab- 
sorbed by dust is re-processed and re-emitted longwards of few 
jum. The Arp220 template is then rescaled in such a way that 
its bolometric emission between 8 and 1000 yum coincides with 
the dust reprocessed luminosity. Finally, the properly scaled 
Arp220 SED is K-corrected and convolved with the MIPS filter 
response to compare with the observed 24 yum flux or flux upper 
limit (see also Berta et al. 2004 ). With this procedure, we have 
been able to use the measured 24 jum flux as a constraint on the 
maximum Av allowed to each photometric redshift solution. 

If compared to other typical LIRG galaxies (e.g. the pro- 
totype local starburst M82), the ULIRG Arp 220 is charac- 
terized by a much more peaked FIR emission and dust self- 
absorption. The choice of such an extreme galaxy as an IR 
template then implies much higher bolometric infrared lumi- 
nosities for a given 24 //m flux, providing an indication of the 
maximun value of Ay consistent with the observed 24 yum flux. 

The results of the SED fitting analysis are reported in 
Figure |8] where for each sources we show on the left panel the 
confidence levels derived from the statistics (averaged over 
all the available free parameters) as a function of redshift z and 
Ay, for the BC03 models (the MA05 model provides equivalent 
results in terms of photometric redshift solutions, see Section 
I4.4l i. Dotted, dashed and dot-dashed curves respectively mark 
the 68%, 90%, 95% and 99.9% confidence levels of the x"^ 
statistics. For each sources, the red shaded area indicates a re- 
gion of the parameter space which is disfavored by our compar- 
ison with the 24 yum flux constraint: spectral solutions falling in 
this range, all heavily dust extinguished, tend to produce 24 yum 
fluxes in excess of the observed values. Although somewhat 
model-dependent and unable to associate a formal significance 
value, our analysis is made relatively robust by our reference to 
the most extinguished object - Arp220 - known at the present 
cosmic time. In particular, we have verified that the adoption 
of other IR spectral templates, like that of M82 or other non 
standard starburt galaxies, would imply somewhat wider ex- 
tension of such "disfavoured" regions (in the sense that even 
lower values of the Av parameter would be inconsistent with 
the measured 24um flux or the adopted upper limits). Note that 
the squared shape of the shaded area is shown here for illustra- 
tive purposes only. 

4.3. Results of the SED fitting 

On the right panel of Fig. [8] the value of the x^ is plotted as a 
function of redshift. Different curves show the result of using 
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different extinction ranges in the SED fitting procedure with 
Hyperz (solid curve: Ay < 6, red dashed curve: Ay < 10). 

As a first result, this analysis indicates the existence of 
multi-modal solutions for the majority of our sample sources 
(14 out of 20, #1, #2, #3, #4, #5, #7, #10, #12, #13, #14, #15, 
#16,#17, #19), even accounting for the 24 yum constraint: gen- 
erally these objects have statistically degenerate fits with either 
a lower-redshift (z < 4) and high-extinguished {Ay ~ 3) SED 
or with a higher-redshift (z > 4) and lower extinction. Sources 
#4 and #19 have best-fit solutions at z >~ 4. 

Six sources (#6,#8,#9,#1 1,#18 and #20) show a single pref- 
erential best-fit solution if we consider the constraint from the 
MIPS emission and if we limit the analysis to the 95% confi- 
dence level. This number reduce to three sources (#6, #8 and 
#11) if this constraint is not taken into account. With the ex- 
ception of the dusty source #9, all these objects favour a low- 
extinction solution: four of them lie at z ~ 4. Only source #6 
favours an higher redshift best-fit (z ~ 8, see dedicated discus- 
sion in Sect. 16.1b . In any case, it should be stressed that, even 
in the case of single best-fit SEDs, the the z ~ 4 solutions are 
only "formally preferred", and alternative lower-z fits cannot 
be ruled-out with any confidence. 

Stellar masses are computed using an adapted version of 
the Hyperz code which performs SED fitting at a fixed red- 
shift. We report in Table [3] a summary of the best-fit parameters 
computed by Hyperz- Only the solutions matching the 24 fim 
constraint are taken into account in the following. The upper 
panel of the table reports the objects with degenerate bimodal 
redshift solution. In this case, for each parameter we present 
the lower redshift primary solution (sol. I) and the secondary 
solution (sol. II) at higher redshift. Moreover, for each physi- 
cal parameter we report a range of values: the value reported 
on the left side of the interval corresponds to the BC03 model 
expectation, while the value reported on the right side of the 
interval corresponds to that of the MA05 model. In the lower 
panel of the table we report the same information for the six 
sources that we consider to have a single photometric redshift. 

For sources with a bimodal behaviour, we compare in 
Figure |9] their observed SEDs (open diamonds) with the cor- 
responding low-z best-fit solutions (left panel, solid red lines) 
reported in Table 3 (sol. I), for the BC03 model. Similarly, in 
the right panel the observed SEDs are shown together with the 
corresponding best-fit high-z solutions (dashed red lines, sol. II 
in Table 3). Similarly, in Figure[TO]we report the best-fit solu- 
tion for the six non-degenerate objects (solutions presented in 
the lower panel of Table 3). 

4.4. Comparison of the Maraston and Bruzual & 
Chariot Models 

As mentioned in the previous Section, we have considered 
two independent evolutionary population synthesis models, by 
MA05 and BC03, in order to check the stability of our results. 
In Table 3 we report the results for both models. The MA05 
model accounts in detail for the contribution of the thermally- 
pulsating Asymptotic Giant Branch (TP-AGB) phase of stellar 
evolution. The TP-AGB phase is calibrated with local stellar 



populations and is the dominant source of bolometric and near- 
IR energy for stellar populations in the age range 0.2 to 2 Gyr. 
The two models have an underlying different treatment of con- 
vective overshooting and Red Giant Branch stars. 

In a recent paper, Maraston et al. (2006) have tested their 
code on a sample of high-redshift early-type galaxies (1.4 < 
z < 2.5), for which the MA05 models indicate younger ages 
(by a factor up to 6) and lower stellar masses (by ~6Q% on 
average) with respect to results obtained with the models of 
BC03. 

For our sample we found that the photometric redshifts ob- 
tained with the BC03 models are very similar to those obtained 
with the MA05 model. The bimodal distributions in the (z. Ay) 
parameter space and the distributions are essentialy equiva- 
lent for the two independent models. A similar situation applies 
for the remaining sources. As can be seen in the second and in 
the third columns of Table 3 (z-phot: sol. I and sol. II), the scat- 
ter in the photometric redshift estimates with MA05 and BC03 
is very small at both low and high redshifts. When considering 
the solution at lower redshift, the mean ratio of the photometric 
redshifts derived from BC03 and MA05, respectively, is 1.04 
with an RMS of 0. 12, and for the solution at higher redshift the 
mean ratio of the photometric redshifts is 1.07 with an RMS of 
0.18. A slightly larger discrepancy is found for the reddening, 
with no clear trends. 

The most interesting physical parameters are the age and 
the stellar mass predicted by the two models. In the lower red- 
shift solution the MA05 model predicts older ages than BC03 
for only 2 out of our 14 objects. For the higher-z solutions 
(including the six non-degenerate objects), 11 out of the 20 
sources with a high-z prediction are older for the MA05. As 
for the stellar masses, we confirm the results of Maraston et al. 
(2006): in the low-z solutions (mean z ~ 3.2) the BC03 models 
predicts in general 50% higher masses on average, and a simi- 
lar trend is found for the high-z (mean z ~ 5.5) fits (~ 60%). 

4.5. Nature of the 24^m Emission 

As discussed in Sect. 13.21 the excess 24 fj.m emission often ob- 
served in the SEDs of our galaxies (Figs.l9]and[T0]i requires the 
presence of ongoing star formation or AGN activity. 

In the extinguished low -redshift case (1 < z < 3), the rest- 
frame MIPS band samples the emissions by warm dust and 
PAH molecules. Both a starburst and an AGN emissions are 
then possible explanations of the 24 pm flux. We decided in 
this case to associate the whole mid-IR flux to a starburst pro- 
cess and adopted the Arp220 SED template (from the spectral 
library by Polletta et al., 2006). As described in Sect.[3T2l this 
template has been used to compute the IR luminosity (8-1000 
p) of each object matching the observed 24 pm flux (dashed 
cyan lines in Figure|9l). 

At higher redshifts (z > 3.5), the low values of the ex- 
tinction favour a scenario dominated by dust-free sources. In 
this case the 24 pm flux has been independently modeled 
with the emission of a dusty torus reproducing the contribu- 
tion of a type-2 AGN (Fritz et al. 2006). Again, in Figures 
|9] and [TO] we show the best-fit solutions at high-z, superim- 



18 



G. Rodighiero, A. Cimatti, A. Franceschini, M. Brusa, J. Fritz and M. Bolzonella: High redsiiift massive galaxies 



posed with the type-2 spectral templates that better reproduce 
the S (lAfim) / S i^fiin) flux ratio (green dotted lines). 
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Fig. 1 1 Distribution of background counts in the full (0.5-7 
keV) band. The histogram shows the distribution of back- 
ground counts in the source cells derived from the Monte Carlo 
method (10,000 trials). The vertical dashed lines show the 3cr 
fluctuation limit: stacked detection with a number of counts 
larger than these values are considered as "real" detections at 
the 99.9% confidence level. The big (red) asterisk represents 
the counts actually detected from the sample of 17 individually 
X-ray undetected sources in our sample. 



5. X-RAY EMISSION AND AGN ACTIVITY 

The best independent evidence of nuclear AGN activity in the 
nuclei of our sources arises if the source is identified as a lu- 
minous (Lx > 10**^ erg/s) X-ray source. X-ray data for a total 
exposure time of ~1 Msec obtained with the Chandra observa- 
tory in the GOODS/CDFS field have been published and made 
publicly available by Giacconi et al. (2002, see also Alexander 
et al. 2003). We have searched for individual X-ray emission 
from the 20 objects in the present sample by cross-correlating 
the IRAC positions of our targets with the positions of X-ray 
sources as catalogued by Giacconi et al. (2002) and Alexander 
et al. (2003). 

We found that three sources (objects #6, #1 1 and #13) were 
individually detected in the X-rays; for all of them the X- 
ray/IRAC positional difference was smaller than 2". The three 
sources share similar X-ray properties: all of them have soft 
X-ray spectra and have been detected only in the soft (0.5-2 
keV) or full (0.5-8 keV) X-ray band, with -13-30 counts, i.e. 
at the limiting flux of the X-ray observation (~ 1 x 10"'^ erg 
cm"^ s '). The large values observed for the X-ray to optical 
flux ratio (FxfPopi) are about 10-100 times higher than gener- 
ally found for AGNs and make these objects "extreme" with 
respect to the overall X-ray source population. The basic X-ray 
properties, as drawn from Alexander et al. (2003), are reported 
in Table |2] 



The low number of counts prevented us to perform detailed 
analyses of the X-ray spectral properties of the sources, espe- 
cially about the amount of gas absorption, given that at z ~ 4 
(~ 8), the 0.5-2 keV observed band corresponds to an intrinsic 
2.5-10 keV band (5-20 keV band) and is thus sensitive only to 
column densities larger than 5x10^^ cm"^ (5x10^-' cm"^). With 
a conservative approach, we derived the intrinsic X-ray lumi- 
nosities, assuming an unabsorbed power-law with F = 1 .8 at 
the best-fit source redshifts, and using XSPEC (Version 1 1.3.1) 
to translate the observed count rates into rest frame 0.5-10 
keV luminosities. The results are reported in Table 2. In all 
the cases, the X-ray luminosities are in excess of lO^-' erg s"', 
therefore suggesting that nuclear activity is fuelling the central 
engine. 

Seven sources with a solid X-ray detection (> 25 counts) 
and undetected at the limit of the GOODS/CDFS observations 
were already reported by Koekemoer et al. (2004, hereinafter 
K04), who named these objects Extreme X-ray to Optical ob- 
jects (EXOs). All of them were subsequently detected in the 
IRAC channels (Koekemoer et al. 2005) with best-fit photo- 
metric solutions in the range 2 < z < 5 (only one being a 
candidate z > 7 AGN). Only one of the three X-ray detected 
objects in our sample, #13, is present in Table 1 of K04. The 
other two escaped the K04 selection because of the low num- 
ber of X-ray counts (< 20). Conversely, of the remaining six 
objects in the K04 sample, five escaped our selection criteria in 
the IRAC 3.6 //m flux (K04 source #4, fainter than 1.8 fiiy) or 
in the /T-band flux (K04 sources #2,3,5,7 have Kab <= 23.5 
in our photometry), while K04 source #1 was not considered 
because of a blending problem in the IRAC 3.6 fim image. The 
fainter AT-band and X-ray fluxes of the X-ray detected objects 
in our sample with respect to the published EXOs fit well in the 
narrow ^T-band to X-ray flux correlation shown in Brusa et al. 
(2005), and can be an effect of a higher redshift origin. 

In order to constrain the average X-ray properties of the re- 
maining 17 individually undetected sources, we have applied 
the "stacking technique", following Nandra et al. (2002) and 
Brusa et al. (2002). For the photometry, a circular aperture 
with a radius of 2" centered at the positions of our sources 
was adopted. The counts were stacked in the standard soft, hard 
and fuU bands (0.5-2 keV, 2-7 keV, and 0.5-7 keV). Extensive 
Monte Carlo simulations (up to 10,000 trials) have been carried 
out by shuffling 17 random positions and using the same pho- 
tometry aperture (2 arcsec). The random positions were cho- 
sen to lie in "local background regions" to reproduce the actual 
background as close as possible. The resulting distributions for 
the trials are shown in Fig.fTTI 

No signal has been detected in the hard band, for a to- 
tal effective exposure time of 14.6 Ms. An excess of counts 
above 3cr of the expected background level was instead clearly 
detected in the soft and full bands. Assuming an unobscured 
F = 1.8 power-law spectrum, the stacked count rate in the 
0.5-7 keV band (2.46e-6 cts s"', or 2.3 10"^^ erg/cm^s) cor- 
responds to an average 0.5-10 keV rest-frame luminosity of 
~ 1.2 x 10"*^ erg s"\ at a z= 4 representative redshift. We have 
verified, by splitting the analysis into various subsamples, that 
this stacked signal was not due to a few objects brighter than 
the average, but was uniformly spread in the sample. 
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Fig. 12 Plot of the ratio of the 3.6 /im to total X-ray fluxes against the 4.5 to /T-band monochromatic fluxes. Our sample sources 
are the red asterisks (the 3 X-ray detected objects) and the open circles (the 17 individually X-ray undetected sources). For the 
latter we report the X-ray to 3.6 ;um flux ratio assuming for all 17 objects the 2.3 10"'^ erg/cm^s total average flux detected 
through the '"stacking"' analysis. These data are compared with data on an X-ray selected sample by Franceschini et al. (2005), 
including type-1 and type-2 AGNs (filled and open squares) and type-2 objects dominated by the host galaxy emission (3-legged 
stars). The lines are the predicted broadband colours as a function of z; the colour for z = corresponds to the point marked 
with 0, and lines are drawn from here to points corresponding to z = 0.5, 1, 1.5, 2, 2.5, 3, and 3.5. From right to left, the lines 
correspond to type-1 quasar and Seyfert-1 (the two red solid lines), Seyfert-2 (the green short-dashed line), and Sb spiral/starburst 
(red long-dashed line). 



We plot in Figure [12] the ratio of the 3.6 yum to total X-ray 
fluxes against the 4.5 to /T-band monochromatic fluxes. Our 
sample sources are the 3 X-ray detected objects (red asterisks) 
and the 17 individually X-ray undetected sources (blue open 
circles). For the latter we report the X-ray to 3.6 yum flux ra- 
tio assuming for all 17 objects the 2.3 10"'^ erg/cm^/s total 
average flux detected through the '"stacking"' analysis. These 
data are compared with data on an X-ray selected sample by 
Franceschini et al. (2005), including type-1 and type-2 AGNs 
(filled and open squares) and type-2 objects dominated by the 
host galaxy emission (3-legged stars). The lines are the pre- 
dicted broadband colours as a function of z for a type-1 quasar 
and a Seyfert-1 (the two red solid lines), a Seyfert-2 (the short- 



dashed line), and Sb spiral/starburst galaxy (red long-dashed 
fine). 



Fig. fT2l confirms that 3 of our high-z galaxies almost cer- 
tainly contain an obscured AGN, because they fall in the AGN- 
dominated region, and in view of their large inferred X-ray lu- 
minosities. The remaining individually undetected objects fall 
in a region intermediate between that of star-forming galaxies 
and of AGNs. Fig.[T2]indicates some probable AGN contribu- 
tions in the bulk of our high-z galaxy sample. These results are 
then not inconsistent with those of Figs. |6] and [T] 
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Table 2 X-ray detected sources in our sample 



# 


X-ray ID" 


A{X - IRAC) 


Counts" 


Counts" 


Counts" 


r b 


best-fit Zphot 


11-4569 
13- 5021 
6- 10945 


217 

205'- 

232"^ 


1.90" 
0.21" 
0.06" 


17. r^-^ 

1 7 4+8:o 


< 15.2 

29 6+^ - 
13 7+'' - 


< 14.2 

< 15.5 

< 12.8 


1.1!;;^ X 10*-' 

2.7+0-5 X io« 
3.6+1-^ X 10« 


4 
4 

8t 



"X-ray data from Alexander et al. (2003) 

''Intrinsic X-ray luminosity at the best-fit photometric redshift, from Table 2 

"^These sources are also present in the Giacconi et al. (2002) catalogue, as XID # 217 and XID #557, respectively, with very similar X-ray 
parameters 

* Following Zheng et al. (2004), Mainieri et al. (2005) associate the observed X-ray emission of #6 (XID #557 in their papers) to a faint optical 
source at a larger distance from the X-ray centroid (~ 1.4"), and derive a photometric redshift z = 1.81. However, the detection of a significant 
2.2 yum signal in correspondance of the IRAC centroid make object #6 the most likely counterpart to the X-ray sources (see Figure [T3]l. 




Fig. 13 Identification of source #6. Left panel: the cutout shows a 10"xlO" negative map of the HST z-band. The black cross 
marks the centroid of the 3.6 yum detected source, while the corresponding IRAC contours are shown in blue (2, 3, 4, 5, 6, 7, 10, 
12, 15cr). The green circle (2" radius) represents the position of the X-ray source as reported by Alexander et al. (AID #232). The 
red circles (1" diameter) indicate the optical sources detected in the field. The thicker red source is the optical counterpart for the 
X-ray source reported by Zheng et al. (2004) and Mainieri et al. (2005). Right panel: a 10"x7" map of the ISAAC /T-band, with 
overlaid the circle indicating the position of the X-ray source (2" radius). There is a significant 2.2 fim signal in correspondance 
with the IRAC centroid, excluding the identification with other sources in the optical map of the left panel. 



6. CANDIDATE MASSIVE GALAXIES AT VERY 
HIGH REDSHIFTS 

6.1. ID-6: A Very Massive Candidate Galaxy atz -87 

We mentioned in Section 14.11 that the SED of source #6 for- 
mally favours a photometric redshift at z ~ 8, with a mass of 
the order of M ~ 8 - lOx IO^Mq. Given the potential relevance 
of this finding, we further analyse this object here. 

Figure [13] details the multi- wavelength identification of the 
source. The Chandra X-ray object associated to this source, 
in particular, has been differently identified by various authors. 
The left panel in the figure shows a 10"xlO" map of the HST 
z-band overplotted to the IRAC contours shown in blue. The 



black cross marks the centroid of the 3.6 fim detected source. 
The (2" radius) double-lined circle represents the position of 
the X-ray source as reported by Alexander et al. (AID #232), 
while the red circles (1" diameter) indicate the optical sources 
detected in the field. The thicker circle marks the optical coun- 
terpart for the X-ray source reported by Zheng et al. (2004) 
and Mainieri et al. (2005). We report in the right panel of Fig. 
[T3] a postage stamp of the /T-band image, showing that there 
is a faint but significant 2.2 /im source in correspondance with 
the IRAC centroid and excluding the identification with other 
sources in the optical map (all at > 1.8 arcsec distances). In ad- 
dition, the X-ray and IRAC centroids are spatially coincident 
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Fig. 14 Two main possible solutions for object ID#6. In the lower panel we consider the primary solution at z ~ 7.7, while in 
the upper panel we show for comparison a secondary solution at z ~ 3.0. The stellar component of the source (solid black lines) 
in both cases correspond to the solution already explored in Section 143] and in Figures l9lfT0l We compare the observations with 
various prototype templates: a type-1 QSO (three dots-dashed green lines, Elvis et al., 1994), the classic type-2 QSO NGC6240 
(dotted red lines, Hasinger et al. 2001), and a typical Seyfert 2 (dot-dashed blue line). The models have been redshifted at the 
corresponding photometric redshifts and normalized to fit the observed 24 yum flux. We have also compared the two solutions 
with the recently discovered most luminous Compton-thick AGN at z ~ 2.5 (SWIRE J104409.95-H585224.8, Pofletta et al. 2006). 
The SED of this object (redshifted and normalized to the 24 fj.m flux) is represented by the filled black circles. 



and clearly indicate the presence of an IR and X-ray source at 
this position which is completely absent in the optical images. 

We have then further investigated the SED properties of this 
object by combining its IR and X-ray information. The com- 
plete observed SED of the source is reported in Figure [141 In 
the lower panel we consider the primary solution at z ~ 7.7, 
while in the upper panel we show for comparison a secondary 



solution at z ~ 3.0. This latter was not formally excluded 
by our analysis, on the basis of the statistics discussed in 
Section l4n a secondary minimum is present at this lower red- 
shift in both panels of Fig. |8] We discuss also this solution in 
the present section. 

In Fig. [14] we compare the SED data for ID#6 with vari- 
ous redshifted prototype spectral templates: a type-1 QSO from 
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Elvis et al. (1994), the type-2 QSO NGC6240 (Hasinger et al. 
2001), and a typical Seyfert 2 galaxy, all spanning the entire X- 
ray-to-infrared wavelength range. For both redshift solutions, 
the models have been normalized to fit the observed 24 fim 
flux. The black continuous-line spectrum from the near-IR to 
the UV corresponds to the integrated stellar component solu- 
tion for the galaxy already found in Sect l4.1l and Figs.l9lfT0l It 
is clear that a purely stellar spectrum cannot explain the large 
observed flux excess at 24yum. In the low-z solution, this ex- 
cess can be explained as due to dust re-radiation from either a 
starburst or an AGN, whereas in the high-z hypothesis the only 
explanation would be a dusty AGN emission. In either case, a 
type-1 QSO spectrum is entirely inconsistent with the data. 

This simple analysis provides interesting hints on the nature 
of source #6: in particular, if we refer to the spectral template 
of NGC 6240 normalized at 24 yum, the observed X-ray flux 
is best explained within the higher-z solution, with a redshift 
of z ~ 7.7. The low-z, z ~ 3.0, fit shows instead a difficulty 
in explaining the X-ray data for both the NGC 6240 and the 
Seyfert-2 templates, since in both cases there is a mismatch of 
the 24//m to X-ray flux ratio, and excess emission would be 
expected in the optical/near-IR. 

A very interesting comparison of the two redshift so- 
lutions is also possible with data on the most luminous 
Compton-thick quasar recently discovered at z ~ 2.5 (SWIRE 
J104409.95+585224.8) by Polletta et al. (2006). The SED of 
this object (redshifted and normalized to the observed 24 jjm 
flux) is shown as black filled circles in both panels of FigfT4l 
Again, for the low-z (lower panel) case, the expected quasar X- 
ray flux would be too faint and undetectable in X-rays. When 
redshifted to z ~ 7.7, instead, this quasar spectrum turns out 
to fit remarkably well the far-lR and X-ray data of source #6: 
the very high redshift brings the unabsorbed hard X-ray spec- 
trum into the observational soft X-ray wavebands. The optical 
and near-lR emissions would in any case be due to the stellar 
component of the host galaxy (solid line, see Sect l4.1l ). 

The analysis presented in this Sect, provides only a quali- 
tative support to the case for a very high-redshift solution for 
source #6, which would imply the existence of an extremely 
massive galaxy (M > 8 x IO'^Mq) at z ~ 8. Clearly, given 
the potentialy extraordinary nature of this object, the lower- 
redshift z ~ 3 solution apparent in Fig. 8 must stifl be consid- 
ered the more conservative conclusion. Unfortunately, proving 
this alternative might turn out impossible before the advent of 
JWST and ALMA. The implications of such a finding will be 
discussed in the next Sections. 

6.2. ID-5: The Mobasher et al. Candidate High-z 
Galaxy 

Our object #5 has been originally identified and analysed by 
Mobasher et al. (2005, object HUDF-JD2), and subsequently 
discussed by Dunlop et al. (2006) and Fontana et al. (2006), 
among others. The availability of extremely deep photometric 
data from UDF and near-IR spectroscopy prompted us to per- 
form a more detailed analysis of this source. To this end we 
adopted a model by Fritz et al. (2007, in preparation; see also 



Berta et al. 2005 and Poggianti et al. 2000). The model spec- 
trum is obtained by summing SSP spectra that are weighted by 
different mass values. Each SSP is extinguished with extinction 
values that are allowed to vary as a function of age (selective 
dust extinction). Dust is assumed to be distributed in a uni- 
form screen (Ry = 3.1). The model computes also the equiva- 
lent widths of all relevant interstellar emission lines. Nine SSP 
spectra have been assumed with ages that were allowed to vary 
from 10*' to ~ 3 ■ 10^, that of the older stars being chosen so as 
to be consistent with the age of the universe. The best fit model 
was found with a analysis after convolution with the filter 
response functions. Furthermore, a constraint was used on the 
Ha emission based on the near-IR spectroscopic observations 
by Mobasher et al. 

Making use of the photometric data by Mobasher et al. we 
explored both the low- and high-z solutions for source #5, run- 
ning our photometric code in the redshift ranges 1 < z < 3 and 
5 < z < 8. In the low-z case, we used as a further constraint the 
Ha line flux upper limit, which we have taken to be 9 x lO'"' 
erg/s from Mobasher et al. 

Our code found two best - fits with photometric redshifts 
z ~ 2.82 and z ~ 6.5 in the two considered intervals. These 
values are fully consistent with our preliminary investigations 
(Sect l4.ll l. The two solutions are shown in Figure [15] The up- 
per panel corresponds to the z=2.82 solution, the lower panel 
to that at z=6.5. The;\f^ statistics sligthly favours the z ~ 2.82 
solution (x^ = 1-27), but the z ~ 6.5 scenario is equally accept- 
able (x^ - 1.82). We have found, in particular, that the lack of 
any Ha emission signal may still be brought into consistency 
with this low-z interpretation if selective dust extinction (i.e. 
extinction as a decreasing function of the age of the various 
contributing stellar populations) is taken into account. 

The z=6.5 fit corresponds to a single-burst solution, with an 
age of 0.3 Gyr, that we already found (Sect. 1431 ). In Figure [16] 
we show the SF history retrieved by the model for the z=2.82 
solution of Fig. [15] In this case the spectrum is reproduced by 
SSP ranging from 2 10^ to 10^ yrs and extinction as a strong 
function of the SSP age. The bulk of the stellar mass would be 
produced between 10^ and 10^ yrs before the observation. In 
view of all these considerations, the lower redshift solution for 
galaxy #5 appears as a more likely explanation. 

7. MAIN IMPLICATIONS 

The results of the analysis presented in this paper are relevant 
in the general framework of galaxy formation and evolution, 
and their main implications are hereby outlined. 

7.1. Massive Galaxies atz ~4? 

Four sources in our sample (#8, #1 1, #18 and #20) show a best- 
fit photometric redshifts around z ~ 4. Solutions at lower red- 
shift produce significantly worse although they are still for- 
mally acceptable. In addition, two more objects (#4 and #19) 
from the bimodal sample present all the acceptable solutions 
at high redshift, with fits of comparable quality at z ~ 4 and 
z ~ 6, and very low extinction. We note here that only source 
#4 is detected at 24/im. 
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Fig. 15 SED analysis for source #5 (Mobasher's et al. source) based on a detailed photometric model. The upper panel corresponds 
to the z=2.82 solution, the lower panel to the z=6.5 solution. The observed photometric data are marked as red-triangles. The 
best - fit model is represented with the blue lines. The convolution of the model in the various photometric bands is shown as 
open blue circles. Note that in this case we used Mobasher's photometry based on the UDF, in order to directly compare their 
results with our modellistic predictions. 
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ysis, which has detected significant average signal (at better 
than 3cr). 

7.2. Extremely dusty starbursts atl < z <3? 
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Fig. 16 Upper panel: Star formation history of source #5 for 
the z - 2.82 solution. The extinction of each SSP populations 
is also reported (green points) on the right vertical axis. Lower 
panel: cumulative distribution of the stellar mass assembled in 
the galaxy as a function of the SSP age. 



All these z ~ 4 candidates share similar physical properties 
with each other. Assuming that z ~ 4 is the correct interpreta- 
tion, in most cases the spectral best-fit solutions correspond to 
either the single-burst model (see Sect. 14.1b or to an exponen- 
tially declining SFR with short decay time-scale t - 0.1 Gyr. 
The typical ages are of the order of ~ 1 Gyr or younger with the 
BC03 models (-1.4 Gyr with the MA05 model). The SEDs are 
then consistent with galaxies at z ~ 4 observed several hundred 
million years after a powerful burst of star formation producing 
stellar masses around 10" Mq with both the MA05 and BC03 
models, and for a Chabrier IMF. 

At lower redshifts, numerous massive (M > lO^M©) and 
old (1-4 Gyr) early-type galaxies have been recently detected 
from A'-band selected surveys in the range 1 < z < 2 (Cimatti et 
al. 2004, McCarthy et al. 2004, Glazebrook et al. 2004, Daddi 
et al. 2005, Saracco et al. 2005, Kriek et al. 2006). Similar de- 
tections at slightly higher redshift have also been recently re- 
ported (z ~ 3.7, Brammer & van Dokkum, 2007) The photo- 
metric modelling generally suggests that these passive systems 
should have formed all their stars at z/ >3 in a short burst. 

Our sources appear to be among the oldest stellar systems 
at z ~ 4, given that the estimated ages are close to the age 
of the Universe at that redshift. The post-starburst nature of 
our z ~ 4 sources and their typical large stellar masses around 
M ~ 1O"M0 suggest that they could be seen as the progenitors 
for the most massive spheroids that are observed around z ~ 2. 
Moreover, as akeady mentioned, the z ~ 4 population is mostly 
undetected at 24 yum, making the passively evolving nature of 
these sources a plausible hypothesis. 

One of this galaxies (#11, see Sect. 5 and FigfT2]i has been 
detected in X-rays and is consistent with hosting an obscured 
AGN. X-ray emissions from all the other 5 objects mentioned 
in this Section have been investigated by us via a stacking anal- 



As mentioned in Sect. 14. 11 the remaining fraction of the sources 
in our sample show a strong bimodality in the photometric red- 
shift solutions. We also showed there that these objects are con- 
sistent with a population of heavily dust-enshrouded starbursts 
at redshift 2 < z < 3, since only few of them are undetected 
at 24 yum . Based on this flux, we have computed the instrinsic 
bolometric luminosity, L/s = L(8 - 1000/im) (see discussion in 
Sect.lO 

The derived luminosities qualify most of these sources as 
ULIRGs (LiR > 10'^ Lo). Assuming that the mid-IR emis- 
sion is mostly contributed by star-formation processes, we can 
translate the IR luminosity into a SFR adopting the standard re- 
lation of Kennicutt et al. (1998). The derived SFRs have a me- 
dian value ~ 3000^^0/3^, even greater than the typical value 
of submillimetre selected galaxies at z ~ 2.2 (~ HOOMo/yr, 
Chapman et al. 2004). The large bolometric luminosities (L/r > 
IO'^Lq) implied by our photometric analysis and the instrinsic 
large extinction required {Ay ~ 2-4), make our sample very 
similar to the optically obscured sources identified by Houck et 
al. (2005) and could represent the most extreme cases of dusty 
galaxies detected until now. 

It seems likely, however, that part of this IR emission may 
be due to an AGN contribution, if we consider that the average 
X-ray signal and the corresponding X-ray-to-3.6 fim flux ratios 
in Fig.[T2]reveal X-ray activity in excess of that expected from 
star-forming galaxies. 

In spite of the intrinsic uncertainties in the photometric so- 
lutions, we have computed the contribution of these sources to 
the SFR density of the Universe by assuming that their far-IR 
emission is mainly contributed by star formation processes. We 
have made us for this of the 1/Vmax estimator (see more details 
in Sect. [8]). Even including all objects with a photometric-z so- 
lution in the redshift range 1 < z < 3.5 and detected at 24 jum 
(#1, #2, #3, #5, #7, #9, #10, #14, and #16), we found that these 
sources missed by optical surveys contributes for only ~20% to 
the global SFR density at 1 < z < 3.5. We obtained a value of 
PsFR-0.026 Ms/yr/Mpc^ (see Perez-Gonzalez et al. 2005 for 
a recent compilation of the Lilly-Madau diagram of IR-selected 
galaxies). 

7.3. Massive Evolved Post-Starburst Galaxies at 

z>6? 

The fourteen sources with bimodal photometric redshift solu- 
tions (upper panel in Table 3) have all a secondary best-fit in 
the redshift range 4.5 < z < 9, see Sect. 14. 11 and Fig. [10] and 
appear to have been detected less than 1 Gyr after a powerful 
burst of star formation producing a stellar mass of the order of 
~ IO^Mq. 

Two sources in particular, (#4 and #19) have formally bet- 
ter solutions at z > 6 and predicted masses M > lO'^M© from 
both evolutionary synthesis models. In addition, source #6 has 
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a best-fit at z ~ 7.7. These three sources have predicted ages for 
their stellar populations of < 500 Myr at their estimated photo- 
metric redshifts of 6.8, 6.2 and 7.7, in the order. Would all these 
sources be at z > 6 galaxies, then these would have formed the 
bulk of their massive stellar populations at very high redshift, 
z > 9, in a period of several tens million years before entering 
a quiescent phase. This would have required an enormous star 
formation activity, with rates of the order of lQ0Q-60Q0MQ/yr. 

Two objects (#4 and #6) have MIPS detections at 24 fim 
and one is bright in X-rays (# 6), which could only be inter- 
preted at such high redshifs as the signature of a luminous type- 
2 obscured quasar within their nucleus (this conclusion being 
proposed also by Mobasher et al. for source #5/HUDF-JD2). 
Mobasher et al. also estimated the size of the dark matter halo 
required to host the stellar mass of source #5/HUDF-JD2 and 
found a value of 4 x IO'^Mq (if the estimated stellar mass of 
the source is 6 X lO^'M© at z ~ 6.5 and a Salpeter IMF is con- 
sidered). 

Assumed such high redshifts and the associated huge en- 
ergetics from both star formation and black-hole gravitational 
accretion, this would have had an important role in the reion- 
ization of large surrouding volumes of the Universe, starting 
the process at redshifts as high as z ~ 15 (Panagia et al. 2005). 
The detection of three comparably massive sources within the 
same small sky region (-130 square arcmin) would support the 
hypothesis that the reionization of the Universe might be dom- 
inated by such massive galaxies. 

However, two of the very high-redshift galaxy candidates 
have lower-z best-fit solutions of comparable quality, which are 
obviously much more likely. Only for source #6 the high-z (z ~ 
8) fit is formally preferred to the low-z one (Sect. lOl i. 

At this stage, and waiting for more decisive future obser- 
vations with JWST and ALMA, we can only set a (stringent) 
limit to the existence of massive evolved galaxies at very high 
redshifts. 

8. CONTRIBUTION TO THE COMOVING STELLAR 
MASS DENSITY AT HIGH Z 

Given the significance of the massive galaxy population at 
z ~ 4 (Sect. 7.1), and even that of the z ~ 6 candidates 
(Sect. 7.3), it is worthwhile to attempt to compare their contri- 
bution to the global comoving stellar mass density with those 
already derived from other independent surveys at those red- 
shifts. We have computed it by using the l/V,„ax estimator (see 
e.g. Franceschini et al. 2006 for an application). Our evalua- 
tion should be considered a lower limit to the comoving mass 
density because our sample is not purely flux-limited (we have 
excluded a priori blended sources and have applied various 
colour limits, see Sect. I2.4l i. 

(A) z ~ 4. We have limited our analysis to the more mas- 
sive galaxies, those with M > lO^M©. For the galaxies de- 
scribed in Sect. 7.1, we have computed for each sources the 
effective co-moving volume of the survey (defined by the sur- 
vey area and the redshift interval, Zm/« < z < Zmax, within which 
each of the z ~ 4 candidates could have been detected, with 
Zmin = 3.7). To determine Zmax, we have taken the best-fit SED 
to each candidates and redshifted it until its flux falls below 
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Fig. 17 The stellar mass function at z ~ 4. The dashed region 
corresponds to the data reported by Drory et al. (2005). We 
have considered the mean value of the mass function in the two 
redshift bins 3 < z < 4 and 4 < z < 5. The lower and up- 
per envelops of the dashed region correspond to the values de- 
rived in the GOODS/CDFS and Fors Deep Field, respectively, 
by Drory et al. (2005). The filled square marks our estimated 
lower limit to the mass function at z ~ 4 for our sample: the 
range spanned by the vertical arrow corresponds to the range 
of values computed with MA05 (lower boundary) and BC03 
(upper boundary) models. 



our limit (S 3.6 - 1 .8//Jy), and taken as Zmax the minimum of the 
corresponding redshift and our upper boundary of the redshift 
interval, z=4.7 (for most of our objects Zmax = 4.7, correspond- 
ing to a survey volume Vmax = 3.9 x 10^ Mpc^ for our adopted 
cosmology). We have derived two estimates of the mass den- 
sity based on the MA05 and BC03 spectral solutions reported 
in Table 3. 

We report in Figure [17] the contribution of our galaxies to 
the stellar mass function at z ~ 4. The shaded region corre- 
sponds to the data reported by Drory et al. (2005), of which we 
have considered an average of their mass functions in the two 
redshift bins 3 < z < 4 and 4 < z < 5. The lower and upper en- 
velops of the region correspond to the mass functions derived 
by Drory et al. in the GOODS/CDFS and FORS Deep Field, re- 
spectively. The filled square marks our estimated contribution 
to the mass function by our z = 4 galaxies. The vertical thick 
errorbar corresponds to our overall uncertainty, in which the 
lower limit corresponds to stellar mass values computed with 
MA05 (based only on 3 objects with favoured photo-z at z=4, 
#11, #18, and #20), while the upper bound is based on BC03 
fits (including 3 sources with favoured photo-z at z=4,#ll, 
#18, and #20, and additional 7 galaxies having bimodal photo- 
z solutions, one of which solutions is within 3.7 < z < 4.7, 
#1,#3,#4, #12,#14,#15 and #19, see Sect|22]above). 

Figure [18] compares our estimated number density at z = 4 
of galaxies more massive than M = 10" M© with literature data 
at different redshifts. The open circle at z = 5.3 is the estimate 
for LBGs with M > 10" Mq recently derived by McLure et 
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Fig. 18 Redshift evolution of the number density of galaxies 
with stellar masses M > 10' 'M©: our data are reported with 
the same symbol as in Figure [17] at z = 4. The open circle at 
Z = 5.3 is the estimate for LBGs with M > lO^M© recendy 
derived by McLure et al. (2006). The filled circles are the esti- 
mated number densities of galaxies with M > 1O"M0 from the 
GOODS/CDFS /T-band selected sample of Caputi et al. (2006). 
We also report the compilation by Drory et al. (2005, open dia- 
monds) in the same mass range. The asteriks marks the number 
density of red and dead galaxies with M > 0.5 x 10"Mo at 
2 < z < 3 detected by Labbe et al. (2005). The solid line shows 
the redshift evolution of the number density of dark matter ha- 
los with masses of M > 2 x IO'^Mq, matching the number 
density of galaxies with M > lO'^M© at z = (adapted from 
Figure 5 of McLure et al. 2006). 



al. (2006), filled circles are from the GOODS/CDFS /T-band 
selected sample of Caputi et al. (2006), and open diamonds 
from Drory et al. (2005) in the same mass range. Altogether, 
our very red galaxies account for a large fraction of the galaxy 
mass density at z = 4 and are among the most massive galaxies 
currently known at such redshifts. 

(B) z ~ 6. As already discussed, our conclusions at higher 
redshift are subject to major uncertainties because of the degen- 
eracy in the solutions for half of our sample (Sections 17. 3l and 
14. 11 1. However, a maximal and minimal case for the contribu- 
tion to the stellar mass density implied by our analysis may give 
interesting insight. For the galaxies described in Section 7.3 
with a photometric solution in the redshift range 5.5 < z < 8.5, 
we then estimated such contribution as for the z ~ 4 case, 
again, we have limited our analysis to the more massive galax- 
ies, those with M > IQ^^Mq. The result of this computation 
is shown in Figure[T8]as the filled square at z ~ 6.7. The ver- 
tical thick errorbar corresponds to the overall uncertainty, in 
which the lower limit corresponds to stellar mass values com- 
puted only with source #6 (the unique object in our sample 
with a clearly favoured solution at high redshift), while the up- 
per bound includes six sources (#4, #5, #6, #10, #12, #16, and 
#19). 
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Fig. 19 Plot of the stellar mass versus photometric redshift, 
summarizing the results of our selection for very high-z galax- 
ies based on selection at IRAC 3.6 yum in the 130 arcmin^ 
GOODS CDFS field, a K-band flux fainter than K ^ 23.5 AB 
mag, and non-detection in any GOODS optical bands. For 1 1 
galaxies shown in filled squares we show the unique photo-z 
solution, for the remaining 9 galaxies we both report as filled 
triangle the lower-z fit and with the open triangle the higher-z 
one. The continuous line corresponds to the mass of a 1 Gyr- 
old SSP detectable at the survey limit. The selection procedure 
clearly tends to detect massive evolved galaxies around z ~ 4. 

The search for z ~ 4 galaxies has so far been mainly per- 
formed with the traditional Lyman dropout colour selection 
technique (Steidel et al. 1999). Lyman Break Galaxies (LBGs) 
are associated with starburst galaxies at high redshifts, identi- 
fied by the colours of their far ultraviolet spectral energy dis- 
tribution around the 912 A Lyman continuum discontinuity 
(Giavalisco et al. 2002). It has been recently suggested that 
60% of the stellar mass at z ~ 5 is missed by the traditional 
drop-out selection technique (Stark et al. 2006), which exclude 
high-redshift galaxies too red in the rest-frame UV to fall under 
the LBG selection (McLure et al. 2006). 

A fraction of our galaxies might represent a complemen- 
tary sample of mature and red galaxies which would, by defini- 
tion, not be selected by the Ly break or other methods making 
use of the optical-UV rest-frame emission. Our results in Figs. 
[T7] and [18] show that the contributions by this new, previously 
unaccounted, population is definitely non-negligible if not to 
dominate the galaxy mass function at the highest z. These ob- 
jects, all escaping detection by published optically-selected or 
K-band selected samples, would almost double the estimated 
value for the high-end of the galaxy mass function (Drory et al. 
2005; Fontana et al. 2006) at that redshift. 

Would this reassessment of the galaxy mass function still 
be consistent with standard ACDM model expectations at such 
high redshifts? Following McLure et al. (2006) and Dunlop et 
al. (2006), we report in Fig. [18] the redshift evolution of the 
number density of dark matter halos with masses of M > 2 x 
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10 Mq (soli line), matching the number density of galaxies 
with M > 10' 'Mq at z = (adapted from figure 5 of McLure et 
al. 2006). The proposed model corresponds to a ratio between 
halo and stellar mass of ~ 20 (see discussion in McLure et al. 
2006). Looking at Fig.[T8]we can see that the number density of 
our massive sources at z ~ 4 is a factor ~2-8 (with MA05 and 
BC03, respectively) lower than the predicted number of dark 
matter halos. Even if our data is only a lower limit, this implies 
that our estimated number density of massive galaxies at z ~ 4 
is still consistent with the current predictions of ACDM models 
for the hierarchical formation of cosmic structures. 

9. DISCUSSION AND CONCLUSIONS 

Recent observational progress has revealed the existence of 
massive structures (galaxies and galaxy aggregates) at high red- 
shift much more frequently than originally supposed. Indeed, 
the problem for any attempts to predict the origin of galaxies 
from first principles is dealing with the very complex physi- 
cal processes involving baryons (e.g. Somerville et al. 2001, 
Somerville et al. 2004). Observations have also shown an un- 
expected trend for the most massive part of the galaxy mass 
function to be put in place first, and lower mass galaxies 
to keep forming stars at lower redshifts (Cowie et al. 1996; 
Franceschini et al. 1998, 2006; Bundy et al. 2005). 

At the present stage, however, the timing for and the phys- 
ical processes accompanying the emergence of massive galax- 
ies with cosmic time, which would be so informative for our 
understanding of galaxy formation, are essentially unknown. A 
problem here is that the most efficient tool to identify very high- 
redshift galaxies, the Ly-dwpout technique, is not sensitive to 
galaxy mass but rather to UV flux. 

The advent of sensitive imagers in the near-IR atmospheric 
JHK bands and particularly of the Spitzer IRAC space facil- 
ity have started to provide new powerful selection tools more 
sensitive to the host stellar mass. Several reports have been re- 
cently appeared about searches for massive galaxies at z > 4 
and the evolution of the stellar mass function (Stark & Ellis 
2006; Yan et al. 2006; McLure et al. 2006; Grazian et al. 2006; 
Fontana et al. 2006). Most of these have in any case exploited 
the Ly-dwpout approach and pushed it to the reddest optical 
bands for the highest-z characterization and used IRAC data to 
constrain the stellar mass. 

Dunlop et al. (2006; see also Mobasher et al. 2005) fol- 
lowed a complementary method of carrying out an extensive 
SED-fitting analysis on large flux-limited samples without pre- 
conceived assumptions about the rest-frame spectrum of the 
candidates, and by using as reference a K-band selection with 
Ks < 23.5. They conclude in favour of a lack of evidence for 
very massive galaxies to be in place at z > 4. 

We follow a similar approach here, but extend it to a galaxy 
selection based on the most sensitive Spitzer 3.6 yum IRAC im- 
ages in the 130 arcmin^ GOODS CDFS, a /^-selection com- 
plementary to that of Dunlop et al. {K > 23.5 AB mag), and 
non-detection in any GOODS optical bands. By these means 
we were aiming at detecting the most massive and highest red- 
shift galaxies in the field to the 3.6 fim limiting flux, and keep- 



ing equally sensitive to dusty star-forming and massive evolved 
galaxies at high z. 

Our results in terms of the stellar mass and redshift for such 
extreme 3.6 /urn selected galaxies are summarized in Figure[T9l 
where we plot the stellar mass versus photometric reshift esti- 
mated from SED fitting. Data for the 6 galaxies with favoured 
single photo-z solution in our discussion, are shown as filled 
squares, while for the remaining 9 galaxies results of both the 
lower-z fits (filled triangles) and the higher-z ones (open trian- 
gles) are reported. 

Given the uncertainties in the photometric redshift solu- 
tions, our results are consistent with previous reports, and point 
towards very few of the galaxies in the field being found at 
z > 6, with only one such extreme case (galaxy #6) being for- 
mally indicated (though a lower-z fit is still acceptable). Two 
other, similarly massive, galaxies are consistent with z > 6 so- 
lution but do not require it. Remarkably, and similarly to what 
found by Dunlop et al., two of the three candidates have solid 
detections in the Spitzer MIPS 24 jum band. However, differ- 
ently from Dunlop et al., we do not consider this as necessarily 
an argument in favour of a lower-z case, because the presence 
of an obscured AGN could easily explain it. Indeed, we have 
found some evidence for (optically) hidden AGNs in the ma- 
jority of our sample of very red high-redshift galaxies from the 
ultra-deep Chandra X-ray data (Fig.fT2li. which adds to the fre- 
quently observed mid-lR 24 yum excess. We concluded from 
this that there might be room for substantial contribution to re- 
ionization to happen in relation with star-formation and AGN 
activity in massive galaxies. 

One major result of our analysis is about the potential ex- 
istence of a candidate population of massive galaxies detected 
around redshift 4. The majority of our sample galaxies (14 out 
of 20) have a photo-z solution at 3.7 < z < 4.7, and 4 of 
them have best-fit solutions in this redshift interval (not for- 
mally unique, however). Hence, in spite of the small numbers, 
a galaxy population undetected in the optical and extremely 
faint in the /T-band appears to possibly dominate the massive 
end of the galaxy mass function at z = 4. These objects, all 
escaping detection by published optically-selected or K-band 
selected samples, would almost double the estimated value for 
the high-end of the galaxy mass function (Drory et al. 2005; 
Fontana et al. 2006) at that redshift. 

Several of these evolved z ~ 4 galaxies (none of the 4 
with robust z = 4 solution, but 4 of the 7 with '"secondary"' 
solution at such z) display strong excess emission at 24 yum. 
This result is similar to that reported by Dunlop et al. (2006) 
for their selected very high-z population. Due to the large red- 
shifts, this would correspond to rest-frame emission at 4-5 yum, 
hence would be difficult to explain purely as dust reprocessing 
by star-forming regions. Again this result requires rather com- 
mon AGN activity in these high-z evolved galaxies. A support 
to this interpretation comes from the deep X-ray data, reveal- 
ing 2 galaxies with bona-fide z = 4 photometric-z to have clear 
AGN-like emissions (with Lx > 10"*^ erg/s), and the remaining 
objects also showing excess X-ray flux (Fig.[T2]l. although at a 
lower level. Some evidence for probable AGN contributions at 
8 and 24 yum was also found directly in the colour-colour plots 
ofFigs.|S|and|7] 
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This widespread association of very high-z galaxies with 
trace obscured AGN activity might bring an interesting confir- 
mation of the emergent view (e.g. De Lucia et al. 2005; Bower 
et al. 2006; Granato et al. 2004) that AGN feedback could have 
systematically influenced the shaping of the galaxy mass func- 
tion during the epoch of galaxy formation. 
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Table 3 Summary of best-fit parameters from Hyperz: photometric redshift, extinction (Ay), age, mass and IR luminosity. The upper panel of the Table refers to those sources 
showing a bimodality in the photometric redshift solutions (14 objects). For each parameter we present the lower redshift primary solution (sol. I) and the secondary solution 
(sol. 11) at higher redshift. The lower panel refers to sources with a single favoured solution. For each physical parameter we report a range of values: the value reported on the 
left side of the interval corresponds to the prediction of the BC03 model, while the value reported on the right side of the interval corresponds to the prediction of the MA05 
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